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Abstract 

All cells derive nutrition by absorbing some chemical and energy resources from the environment; 
these resources are used by the cells to reproduce the chemicals within them, which in turn leads 
to an increase in their volume. In this study we introduce a protocell model exhibiting catalytic 
reaction dynamics, energy metabolism, and cell growth. Results of extensive simulations of this 
model show the existence of four phases with regard to the rates of both the influx of resources 
and cell growth. These phases include an active phase with high influx and high growth rates, an 
inefficient phase with high influx but low growth rates, a quasi-static phase with low influx and 
low growth rates, and a death phase with negative growth rate. A mean field model well explains 
the transition among these phases as bifurcations. The statistical distribution of the active phase 
is characterized by a power law, and that of the inefficient phase is characterized by a nearly 
equilibrium distribution. We also discuss the relevance of the results of this study to distinct states 
in the existing cells. 

PACS numbers: 87.17.Aa, 89.75.Fb, 82.39.Rt 
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I. INTRODUCTION 



To unveil the minimal requirement of life, we need to first understand how an ensemble 
of a variety of catalytic reactions generates a stable cell that grows to reproduce itself. A 
cell in any form — right from its prototype form at its inception to its currently evolved 
form — receives an influx of some nutrient chemicals; the cell then effectively uses the energy 
obtained from the decomposition of these nutrient chemicals to develop useful catalysts and 
membranes. It should be noted that the chemicals and energy obtained from the nutrients 
through catalytic reactions are used for cell growth and cell replication. 

Given the fact that cell growth requires a continuous flow of energy and nutrients, it can 
be said that cell growth will not be sustained if the cell is in thermal equilibrium. Indeed, 
the relevance of non-equilibrium open systems to life processes has been discussed over 
several decades . On the other hand, in a Carnot cycle, the maximal thermodynamic 
efficiency is achieved under an adiabatic condition, where a weak influx flow is generated 
under a condition close to the equilibrium condition. In this case, a cellular state far from 
equilibrium may not be appropriate for efficient thermodynamic processes. Therefore, it 
is pertinent to determine the condition under which highly efficient cellular growth will be 
achieved f 

This problem regarding efficiency in cell- volume growth is a general problem so that it is 
essential to understand the primitive form of a cell, or a protocell. In this case, one needs 
to understand a universal property of a growing cell following catalytic chemical reaction 
processes upon external nutrients. On the other hand, this issue is also important to the 
present cells. In this case, it is likely that the rate of conversion of chemical resources to 
increase the cell volume depends on the current cellular state, or one may distinguish distinct 
cellular states on the basis of the efficiency of conversion. 

With regard to such cellular states, it is a known fact that a cell undergoes several 
"phases" depending on culture conditions, where each phase is characterized by a distinct 
growth rate [4 (J. The appearance of these phases depends on the nutrient condition; more- 
over, each phase represents a distinct cellular state whose chemical composition differs from 
that of any other state. For example, in the case of bacteria, a high growth rate is achieved 
in the log phase, which is characterized by exponential cell growth under a nutrient-rich con- 
dition; having said that, bacteria also undergo some secondary phases such as a stationary 



phase characterized by suppressed growth and a dormant phase characterized by suppressed 
metabolism. These observations imply the existence of qualitatively distinguishable cellular 
states. With this background, is it now possible to understand the existence of such distinct 
phases, characterized by distinct growth and influx rates, in a cell? The purpose of the 
present study is to introduce a simple "protocell" model having distinct cellular states and 
analyze efficiency for cellular growth depending on each state. 

To investigate the origin of multiple cellular states and underlying energetics theoreti- 
cally, a cell model including metabolism and growth is needed. There are extensive studies 



conducted on protocell models that are based on catalytic reactions |7H17| to determine 
possible conditions that facilitate the reproduction of cells. Previously, Furusawa and one of 
the authors (KK) studied a protocell model consisting of a large number of chemicals that 
mutually catalyze each other, in order to understand the manner in which a cellular state is 



realized 
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15|. Random catalytic networks were used to identify the condition for steady 
cell growth, where we found a universal statistical law on the abundance distribution of 
chemicals. Interestingly, in spite of the simplicity of the model with just catalytic reaction 



networks, this statistical law was also confirmed over most of the existing cells [12 |. 

However, in the model, energetics was not discussed. It was assumed that the energy 
required for catalytic reactions was supplied automatically, and also that the growth rate 
of the cell was equal to the inflow rate of the resources. Furthermore, multiple distinct 
cellular states were not observed with regard to this model. Therefore, the growth rate was 
uniquely determined by the nutrition condition, and the efficiency of cell growth depending 
on cellular states could not be studied at all. 

To study the cell-state dependence of growth efficiency, we extend here the above- 
mentioned protocell model to include energy metabolism. For this purpose, we first need to 
consider reversible reactions, with the reaction rate constants given by the chemical poten- 



tials of molecules to be consistent with thermodynamics 18|, |l9). This consideration is in 
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W22j, 



contrast to that pertaining to most models of catalytic reaction networks 
where only irreversible(uni-directional) reactions are considered. Second, we introduce spe- 
cific chemicals for "energy currency molecules" such as adenosine triphosphate (ATP), to 
derive energy from the reactions or to supply energy to the reactions. This energy currency 
component is assumed to be used for cell growth and is synthesized by other catalytic re- 
actions. In fact, the cell growth usually requires some high-energy chemicals , the synthesis 



of which requires energy currency molecules. For example, the synthesis of DNA requires 
ATP consumption. 



In the case of our model, first we consider that there is an influx of nutrient resources. 
These nutrients are used in the synthesis of a variety of catalysts through reactions given 
by catalytic reaction networks. An energy currency molecule (e.g., ATP) is synthesized 
through such reactions, while the cell volume increases by consuming the energy generated 
by this molecule. Further, the metabolic efficiency of this model, i.e., the rate of cell growth 
for consumption of nutrients, depends on the internal cellular state, i.e., the concentrations 
of the catalysts. Hence we can study the relationship between the influx rate and the cell 
growth rate, depending on the cellular state, and examine the possibility of existence of 
distinct cellular states with unique metabolic efficiencies and cell growth rates. 



The present paper is organized as follows. In the next section we introduce a model 
consisting of a catalytic reaction network with energetics, and cellular growth. In Sect. 
Ill, direct simulation of the model and the mean-field analysis are presented, where we 
report the existence of different phases with distinct cellular growth speed and efficiency, as 
well as different statistical characteristics of chemical abundances. Section IV is devoted to 
summary and discussion. 

4 



II. MODEL 




FIG. 1. Schematic representation of our "protocell" model. It consists of an encapsulated chem- 
ical reaction network within a membrane. Green circles indicate catalytic components and i,j, k 
are the component indices; purple circles indicate energy currency molecules. Arrows across the 
membrane represent material exchange by diffusion through the membrane, whereas arrows within 
the membrane represent chemical reactions. Conversion from energy currency component A to B 
is used for the synthesis of membrane molecules resulting in an increase in the cell volume. 

Here, we consider a very primitive "minimal" protocell, compartmentalized by a mem- 
brane, within which encapsulated catalytic reactions progress (see FigJT]for schematic repre- 
sentation). This catalytic reaction is necessary for the synthesis of all the chemicals within 
the cell, and also brings about cell growth. For this synthesis, some supply of matter and en- 
ergy is generally required. For the former, we consider some nutrient chemical species, from 
which catalysts are synthesized. For the latter, a common molecular species that works as 
an "energy currency", say ATP, is assumed; it is required for many endothermal reactions 
to progress. There are influxes and effluxes of these and other chemicals. As the num- 
ber of molecules increases, the cell volume increases (as a result of synthesis of membrane 
molecules) and, accordingly, the concentrations of all chemicals within the cell decreases. 

For system variables, we consider the following: (i) A set of concentrations of catalytic 
components whose concentrations are given by (xq, x\, ...xn-i)- These chemicals are syn- 
thesized by catalytic reaction(s), and other metabolites are assumed to exist constantly, 
hence they appear as constant parameters in the reaction, (ii) An energy currency molecule. 
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With ATP and ADP in mind, we assume that the conversion of A to B supplies energy (A 
corresponds to ATP, and B to ADP). The concentration variables of these are (jja, Ub)- (hi) 
The volume of the protocell is set at V, and increases or decreases with time. 

Catalytic reaction coupled to the conversion of energy currency molecules: We consider 
a reaction process from % to j catalyzed by /, coupled with the conversion of energy currency 
molecules. This reaction is given by 

i + l + A(B)<*j + l + B(A), (1) 

where i,j,l are reactant, product, and catalyst, respectively. By denoting the stan- 
dard chemical potentials for the catalytic components and energy currency components 
as n®, fi%), the reaction rate constant satisfying the detailed balance condition is 
determined as 

hj = A i:j min[l, e-^+^Vr^-zV))^ (2) 

where /3 is the inverse temperature and A^ = A,-j[19j]. The reaction rate for % to j is then 
represented as kijXiXiyA(B)- Note that the reaction progresses to reduce the free energy of the 
associated chemicals, not the energy. From this point of view, the assumptions of the mass- 
action law and the rate constants correspond to taking the form of the chemical potential as 
fii = /i°+(l//3) logXj. In fact, the net reaction flux from i to j, i.e., x;(/c i jXj|/A(_B)—/cji^ : ,|/s(A))) 
always reduces the free energy per unit volume g = XiHi + Xj\ij + DaHa + Vb^b- 

We consider a variety of reactions. The whole set of reactions constitutes a network given 
by the tensors CA{i,j, and Cs(i,j, I), which represent the structure of a catalytic reaction 
network, defined by 

{1 (If reaction i + l + A(B) -)■ j + 1 + B(A) exists) 
(3) 
(otherwise) 

The nutrient chemical is assumed to diffuse through the membrane, 
and its external concentration is set at X Q , hence the efflux rate is given by D (X — xq), 
with the diffusion coefficient Dq. For simplicity, other components except the chemical N — 1 
are assumed to be impermeable, i.e, Di = 0. The external concentration X^v-i is set much 
smaller than Xq, hence the cell takes as nutrient inflow, and diffuses out N — 1 as waste. 
However, this condition is not essential and can be relaxed. Several other DjS could be set 
at non-zero, and the result would be qualitatively unchanged. 



Cell growth progresses as a result of membrane increases (and the influx of water in 
proportion). Instead of explicitly introducing a membrane molecule, here we simply note that 
membrane growth needs some energy currency consumption. Assuming that the precursor 
molecules for the membrane are abundant in the medium, we consider a phenomenological 
"growth reaction" that progresses accompanied by the reaction A — >■ B. We assume that the 
rate of the growth reaction also satisfies the detailed balance condition. Then, defining P 
and jjp as the concentration and standard chemical potential, respectively, of the membrane 
precursor, and \lm as the (effective) chemical potential of the membrane, we write the rate 
of the growth reaction as 

k{Py A e^ + ^ ] - y B e^ M+ ^) (4) 

where k is a rate constant. In addition, we assume that the concentration of the precursor 
is fixed. Hence, the volume growth is simply given by 

~ = -yk g {v A - yB e*) (5) 

where e (= fi M + fi^ — fi p — (1/ (3) logP — fi° A ) is the energy required for the process, k g 
(= kPe^^ ^^) is its rate constant, and 7 is the conversion rate from the membrane 
molecules to the cell volume. 

Here, P, fip, /zm, fi% /^s, k are direct independent parameters controlling the reactions, 
and kg, e are functions of these parameters. However, the behaviors of the model as dynam- 
ical systems are well represented by choosing k g , e (and X ) as relevant control parameters, 
hence that we use these as the relevant parameters in the present model. 

As a result of the growth in cell volume V, given by eq. ([5]), all the chemicals are diluted 
accordingly. This dilution leads to a decrease in the concentrations of all the chemicals in 
proportion to themselves, i.e., Xi^^r- As a result of this dilution, it is expected that the total 
concentration of the energy currency components (i/a + Db) would decreases, and, without 
a supply, the system would cease to grow. In reality, such molecules would also be diffused 
from the outside, say phosphorate molecules, the material components for ATP and ADP. 
Hence, the energy currency molecules are assumed to flow in with increasing the cell volume 
so that Da + Vb is constant. Consequently, the two concentrations are not independent. We 
therefore only need to consider the time evolution of jja- 

Taken together, the time evolution of the concentrations of the catalysts XjS and energy 
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currency ija are given by 



— ^ ] I ^ ^ C n (j,i,l)xi(kjjXjy n kijXiDm) 

(n,m)=(A,B),(B,A) 

1 

+ A(^i -Xi) - X iy-jf ( 6 ) 



= -&<?(^ - VB^ e ) + ^2C B (i,j,l)xi(k ij x i y B - k^x^A). 

First, all of the coefficients of the chemical reaction rate (A^s) are set at 1 , for simplicity. 
Secondly, the standard chemical potentials of the catalysts (//°s) are sampled from a uniform 
distribution from to 1. Thirdly, if the component is a nutrient, the energy currency A 
is assumed to be synthesized when the component is degraded, by appropriately setting 
the directions between the catalytic reactions and the conversion of the energy currency 
component. For all other reactions, the coupling directions have been determined randomly. 

Finally, we make a cautious comment on the effect of the influx of the energy currency 
molecules on the intracellular energy currency composition. The compositions of the energy 
currency molecules A and B outside the cell are expected to be in chemical equilibrium. 
Then, when the energy currency molecules in such a composition flow in with cell growth, the 
compositions of the intracellular energy currency components will be equilibrated. 

We can treat this effect by (i/a — Va)v^ on ^° the ti me evolution of i/a, where y A q is 
the concentration of i/a at equilibrium, and we carried out simulations for such a model. 
However, we confirmed that the results discussed in the present paper are qualitatively the 
same, independent of this effect; hence we neglect this effect for simplicity. 

III. RESULTS 

A. Emergence of states with different metabolic efficiencies 

We simulated the present model using several parameter values and taking a variety 
of catalytic reaction networks. Each network was generated randomly with a given total 
reaction number G under the constraint that it is percolated. From these simulations, we 
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found common behaviors in the present protocell model, as long as the reaction path number 
was sufficient. Indeed, as long as the average path number per chemical, i.e., K = G/N, 
satisfies K > K* ~ 8, the common behavior is observed, independent of the network. If K 
is smaller than that, the behavior often depends on the network, which we do not discuss 
here. 

For most of these cases, the chemical concentrations and growth rates converge to a fixed 
point after some relaxation time. If the growth rate is positive, a fixed point represents 
a cellular state with steady growth. If it is negative, the cell shrinks, and will ultimately 
collapse. We are interested in the former case here. Since the approach to fixed points is 
quite common, we focus on the properties of such fixed points, without going into details of 
the transient relaxation process. For most parameter values, the stable fixed point is unique, 
but at some parameter values, there are two stable fixed points, leading to bistability in the 
cellular state;this is discussed below. 

In cell state characterization, the growth rate A = and the influx J in = D (X — x ) 
are important order parameters for characterizing the protocell; the "metabolic efficiency" 
77, defined as r\ = \/Ji n , is another basic order parameter. 

With regard to the internal properties of a protocell, the density m is defined as the total 
concentration of catalysts, as m = X^ 2 -*- The density is an important order parameter and 
is related to influx, efflux, and growth rate. In a steady state of the cell, there has to be a 
balance between material influx and dilution by cell volume growth so that the density is 
constant. Hence, we get 

Am = J in - J out (7) 



which gives the relationship among the growth rate, density, influx, and efflux. Then, the 
efficiency r\ is represented by i] = (1 — J out / 1 J in ) /m. The three order parameters (A, Jj„,m) 
are sufficient to specify a cellular state macroscopically. 



9 



1 c i i i r| — i^^^^^ ^B 1000 




(a) (b) (c) 

FIG. 2. Dependences of the steady state properties on growth reaction constant k g for Xq = 0.5 
(blue) and Xq = 10 (red/purple), using a catalytic network with N = 24, K = 10. The parameters 
were set as X N ^ = 0.01, D = 1.0, 7 = 1.0, = 1.0, e = 0.0001, = ^ = 1.0, and y A + VB = 1.0. 
(a)Growth rate:A = p*^-; (b)influx: Jj n = Dq(Xq — xq); and (c)density:m = Yli x i- 

Fig. [2]shows the dependences of the growth rate A, influx J in , and density m on the growth 
reaction rate constant k g for a certain catalytic reaction network with iV = 24, K = 10. 
(Indeed, this behavior is commonly observed, independent of the networks.) When Xq is 
large (Xq = 10), it is found that there are two branches of fixed points with different growth 
rates, as shown in Fig. [2(a). In this case, these two fixed points are bistable when k g takes 
an intermediate value. The transition between the two states with the increase or decrease 
in kg occurs with hysteresis. 

These two states with higher/lower growth rates correspond to lower/higher influx and 
density, respectively, as shown in Figs. E^b) and [2(c). Here, however, the dependence of the 
influx on k g is much weaker, and there is not a large gap between the influx values of the two 
states. Thus, the difference between the growth rates of the two states directly represents 
the difference in metabolic efficiency rj. This significant difference in the efficiency reflects 
the internal state of the cell, as discussed later. 

Fig. [2] also shows that only one state with the intermediate growth rate exists when the 
external nutrient concentration is small (Xq = 0.5). This state has a very low influx value, 
and accordingly, a high efficiency and is therefore distinguishable from the above two states. 
When Xq decreases, we observed that the efficiency increases, which is a common behavior, 
independent of the values of k g . 

Besides the states shown in Fig. [21 another state with a negative growth rate appears 
when k g is large enough. In this case cells shrink, and hence that the state is regarded as a 
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"death" state. We note that the "death" state coexists with another growing state if Xq is 
not small, but, for small Xq values ( say smaller than 0.1), it is the only attractor. 

Summarizing these distinct states, we have plotted the phase diagram against k g and 
Xq as shown in Fig. El Each "phase" here is defined and characterized as in Table [H In 
particular, as expected from the plot in Fig. [2J the range of k g values in which phases I and 
II coexist narrows continuously as Xq decreases, and finally disappears. 
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FIG. 3. Phase diagram of the model on the Xq — k g plane. Each phase (I, II, III, IV) is defined as 
in Table [H Solid lines represent discontinuous transitions between phases; broken lines represent 
continuous transitions. 



Phase 


A 


Jin 


V 


Distribution of abundances 


I: Active phase 


High 


High 


Medium 


Power law 


II: Inefficient phase 


Low 


High 


Low 


Close to equilibrium 


III: Quasi-static phase 


Low 


Low 


High 


Mixed 


IV: Death phase 


Negative 









TABLE I. Definition of phases and characterization by order parameters, and the statistical 
distribution of abundances of chemicals to be discussed in Sec. IIII CI 



B. Mean field analysis 



The existence of each phase and the structures of the transitions between phases are ex- 
plained by carrying out a simple mean field approximation for our model. To construct the 
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mean field equations of our model, we consider the average concentration x of chemicals over 
all chemical components (i = 0, 1, N — 1), i.e., x = m/N, and replace the rate equation by 
this concentration. To get this mean field approximation, the catalytic network is approx- 
imated by a fully connected network with weak paths so that it has the the same average 
path density p in the original model (= K/N). Also, the standard chemical potentials p® 
of all the chemicals are set to be identical (to their mean value), and the standard chemical 
potentials of the energy currency components p° A , p° B are the same. In addition, by assuming 
that diffusion through the membrane is fast enough, the influent component concentration 
x and effluent component concentration x^-i are regarded to be equilibrated with their 
external values, and they are replaced by X and X^-i respectively. 

Combining these approximations, we obtain the following temporal evolution equations, 
which consist of the concentration variables of the catalyst x and the energy currency 
molecule y A , as 

xk g -/(y A - y B e 0t ) (8) 
y B )(x - X A r_ 1 ) 

(9) 

As with the direct numerical results, the existence of two branches of fixed points is 
confirmed for large Xq, with hysteresis as a change of the growth reaction rate constant k g . 
In Fig. HJ^a), the growth rates are plotted against k g for Xq = 0.5 and 10, which are obtained 
from the mean field equations, eqs. ([8]) and (Q. Comparison with the results in Fig. |2]^a) 
indicates that the obtained dependence of the two fixed point values on the parameter k g , 
and the existence of hysteresis for large Xq, agree well qualitatively with the results from 
the original model. 

In addition, the death phase, IV, with a negative growth rate, is also obtained from the 
mean field equations. Again, for a medium value of Xq, this fixed point coexists with that 
for phase I or III. In contrast, with decreasing k g , the region of coexistence narrows and 
finally disappears. 

By summing all the solutions, a phase diagram is obtained from the mean field equations, 
as shown in Fig. H](b);the diagram agrees rather well with that in the original model(Fig. 

H. 
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dx i 

— = px{X y B - xy A ) - -px(y A + y B )(x - X N ^) 



*M = p(N- 2)x(X y B - xy A ) - \p{N - 2)x(y A 
at 2 

- k g (y A - y B e^) - p(N - 2) 2 x\y A - y B ). 




10 



100 



0.01 



10 



0.01 



0.1 



10 



100 



kg 



(a) 



(b) 



FIG. 4. Behaviors of the mean field equations (eqs. ([8]) and ([9])). The other parameters are set to 
be identical with those adopted in Fig. [2J (a) Dependence of growth rate on k g , obtained from the 
steady states of the mean field equations, for X$ = 0.5, 10. Colors represent different Xq values, (b) 
Phase diagram from mean field equations on the Xq — k g plane. Each phase is defined in the text 
as given in the original model. Solid lines represent discontinuous transitions between the phases 
with saddle-node bifurcations, and broken lines represent continuous transitions distinguished by 
changes in the growth rate and influx without bifurcations. 

From the mean field equations, one can see that the mean catalyst concentration x (or 
the density m = Nx) has a negative effect on the growth rate A ; this introduces bistability 
in phases I and II. The third and fourth terms in the right-hand side of eq. (Q indicate the 
consumption of A by the growth reaction, and the degradation of A by catalytic reactions, 
respectively. The dependence of each term on catalyst concentration x shows that the 
degradation of A increases faster than consumption, so that the allocation of energy currency 
to the growth processes decreases with increasing x. Hence, for large x, the abundance of 
A and the growth rate are suppressed. Consequently, the dilution effect, represented by the 
third term in the right-hand side of eq. jSJ), does not work effectively, and x is sustained 
to take a large value. On the other hand, when x is initially small, it continues to take a 
small value, as it is suppressed by dilution. As a result, two states with different growth 
rates coexist. Noting that m = xN in the mean field approximation and replacing x by 
the density m, the above argument is also applied to the original model. To sum up, there 
is mutual inhibition between the growth rate and catalyst density. This mutual inhibition 
leads to bistability. Note, however, that for this mechanism to work, the density m should 
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be sufficiently high, so that bistability is not present for small Xq. In this case, a single 
phase, III, exists against the change in k g . 

Since the mean field equation involves just two variables, the transitions among each 
phase are analyzed as bifurcations of fixed points. From the standard analysis of the flow of 
a two-dimensional state space and the linear stability of the fixed points, it is straightforward 
to show that there are saddle-node bifurcations at the transitions from phase I to II and II 
to I, with kg as a bifurcation parameter. These two bifurcation points correspond to the 
boundaries between the phase I and the bistable region, and the phase II and the bistable 
region. As X decreases, these bifurcation points approach and finally collide and disappear. 
This disappearance corresponds to the transition to phase III. 

Transition to the death phase, IV, has the same bifurcation structure. There are saddle- 
node bifurcations at the transition between the death phase IV and the growing phases with 
X as a bifurcation parameter. These bifurcation points approach and finally collide with 
each other and disappear as k g decreases. 



C. Statistical natures of internal states 



So far, we have shown the distinct phases characterized by "macroscopic" variables, such 
as the growth rate, influx, and total density of chemicals, which are well represented by the 
mean field equations. In a direct simulation involving many chemicals, however, there are 
'microscopic' variables such as the catalyst concentrations (xq, X\, ...xn-i)- In our model, 
there is no specific meaning to each chemical, hence the distribution of abundances is a 
relevant characteristic of a cellular state. For this statistical property of the concentrations, 
we take a rank-ordered concentration distribution, i.e., we plot the abundances in order of 



magnitude, following earlier studies [1 21] . As a result, it is revealed that each phase has a 



quite different statistical property as discussed below. 
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(e) 

FIG. 5. Statistical nature of the model. A network with N = 2000 and K = 20 was used. 

(a) Rank-ordered concentration distributions in phase I for Xq = 0.1,0.3,0.6. Colors represent 
different values of Xq. Here k g = 10.0, and the other parameters are same as those in Fig. [2j 

(b) Distributions in phase II. The network and parameters are the same as those in (a), except 
kg = 0.1. (c)A distribution in phase III. The network and parameters are the same as those in (a), 
except for Xq = 0.02, k g = 5.0. The inset shows a log-log plot of the distribution up to rank 30. 
(d)Dependences of the growth rate on Xq in phases I and II. Colors and shapes of points represent 
data corresponding to the plots in (a) and (b) sharing the same color and shape, (e) Dependences 
of Kullback-Leibler divergence from the distributions to the equilibrium compositions on Xq . The 
inset shows the dependencies on k g . -, p- 



Figs. EI a), (b), and (c) show rank-ordered concentration distributions in phases I, II, 
and III, respectively. The ordinate indicates the concentration of catalysts and the abscissa 
shows the rank determined by the concentration. For phase I, we observed that the slope of 
the rank-ordered concentration distribution increases with increments of X , and reaches a 
power- law distribution with exponent -1. However, we note that this power-law distribution 
with an exponent -1, in other words, Zipf's law, can be achieved only when the catalytic 
reaction network is sufficiently sparse (approximately K < 40 for iV = 2000). For a dense 
network, although the slope of the distribution increases with Xo, the distribution remains 
curved in the log- log plot, and thus does not fall on the power law. The change in the slope 
is accompanied by a change in growth rate, which, for reference, is plotted in Fig. Eld) 
with increasing X for phase I. Each point in the figure corresponds to each distribution in 
Figs. Ela) and Elb) of the same color and shape. With increasing the growth rate, the slope 
of the distribution also increases and approaches the power law with exponent -1, near the 
optimal growth rate. With further increases in Xo, however, the growth rate decreases until 
it reaches zero. In this region before the arrest of growth, the concentrations approximately 
obey a power-law distribution with exponents -1 or slightly higher (data not shown). In 
this region, however, the steady state is often no longer stable, and the growth rate shows 
oscillations with time. The plot with the downward arrow in Fig. Eld) denotes a regime 
with such oscillations. 

Fig. Elb) shows rank-ordered distributions in phase II, for several values of Xo by a 
semi- log plot. In contrast to phase I, the plot is linear on a semi-log scale, and the slope 
is independent of the values of Xo. In other words, the distribution is exponential, which, 
as shown below, comes from the fact that the concentrations have near equilibrium com- 
positions in these states, defined by the standard chemical potentials of the catalysts as 
x eq oc exp(— /3/i°). Because of the uniform distribution of the potentials the rank-ordered 
equilibrium distributions take the form of exponential distributions. 

To confirm that the distributions in phase II are actually close to the equilibrium dis- 
tribution, we measure the distance from equilibrium by calculating the Kullback-Leibler 
divergence from the equilibrium distributions of the concentrations (= J2iPi ^°ziPi/pT)i 



where = xf 1 ' / V, a^ ei ^)|23| as plotted in Fig. Ele), against the change in X ; the 
change of the divergence against k g is shown in the inset of Fig. Ele), ranging between phase 
I and II. It is shown that the divergence in phase II is quite low compared to that in phase I 
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and decreases with increments of Xq, to approach equilibrium. Although the emergence of 
exponential distributions itself depends on the assumption of uniform distribution of /x°s, the 
decreases in the divergence are independent of assumptions on the distribution of /i°s. Thus, 
phase II can be characterized by the distribution close to the thermal equilibrium. As shown 
in the inset of Fig. 15(e), this branch connected continuously from the equilibrium shows a 
discontinuous jump to phase I, where the composition of the concentration distribution is 
far from that at thermal equilibrium. 

Thus, phases I and II are characterized by distributions with power-law distributions and 
close to equilibrium, respectively. On the other hand, the states in phase III are mixed, 
as shown in Fig. |5Jc). For chemicals with smaller abundances the rank-abundance plot 
follows exponential decrease and is near equilibrium, as in phase II. The concentrations 
of abundant chemicals, however, deviate distinctly from the equilibrium distribution. The 
rank-abundance plot for abundant chemicals follows a power law, as in phase I, while the 
fraction of such chemicals that deviates from equilibrium distribution increases with k g . This 
partial deviation from equilibrium distribution characterizes phase III. 

IV. SUMMARY AND DISCUSSION 

We have studied a simple "protocell" model, in which chemicals are encapsulated within 
a membrane. The protocell synthesizes these chemicals by using the nutrients it derives from 
the environment ;the protocell also utilizes the energy obtained from intracellular reactions 
to increase its volume. These reactions are catalyzed by other chemicals that are also syn- 
thesized within the cell, and thus, the total set of reactions constitutes a catalytic network, 
which was randomly chosen in this study. From the results of simulations of a variety of 
networks, we have obtained the steady state in which a cell grows with a constant rate. 

We have carried out simulations by varying the values of the concentration of the nutrient 
outside the cell, i.e., X , and the rate constant of the phenomenological growth reaction, i.e., 
kg. Results of these simulations show that three distinct phases with qualitatively different 
cellular states exhibiting steady growth exist, depending on these two parameters. When 
the external nutrient concentration Xq is high, phase I appears for a large value of /c 9 ;this 
phase is characterized by high cell growth and influx rates. On the other hand, when the 
value of k g is small, phase II appears, characterized by a low growth rate and a high influx 
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rate. At an intermediate value of k g , these two phases coexist. When Xq is low, the influx 
rate decreases, and this regime where phases I and II coexist narrows and then completely 
disappears. After the disappearance of phases I and II, phase III comes into existence;this 
phase is characterized by low growth and influx rates, regardless of the values of k g . Further, 
this phase is characterized by a high efficiency of conversion of nutrients, which enable an 
increase in the cell volume. In addition to these states characterized by a positive growth 
rate, we also found phase IV, termed the "death" phase; this phase is characterized by 
negative growth rates for low Xq. When k g is high, it is found that both the death phase 
and the growth phase I or III coexist. 

These transitions between the phases are explained using the mean field equations of 
the model, which are obtained by assuming that the path density is sufficiently high. The 
mean field equations account well for the transition between phases I and II, as well as those 
between phase I or III and IV as saddle-node bifurcations. The success of the mean field 
equations that are independent of each network also explains why the behaviors observed 
here are independent of the choice of networks. The analysis of the mean field equations 
enables us to understand the bistability of two phases as mutual inhibition between the 
cell volume growth and the concentration of intracellular chemicals. Because this mutual 
inhibition does not depend on the details of the model, we expect that this bistability will 
be a general feature of any cells in which the rate of conversion of nutrient for increasing 
cell volume can vary depending on the composition of chemicals in the cell. 

In the original simulations of many chemical components, we have also examined the 
statistical distribution of the abundance of chemicals. To characterize each cellular phase 
statistically, we adopt the distribution of abundances. In phase I, the abundance of chemicals 
differs from species to species, spreading over several orders, and this distribution is far 
from that expected from the equilibrium distribution defined by chemical potentials. As the 
growth rate increases, the abundance distribution of chemicals approaches a certain power 
law, known as Zipf's law. On the other hand, the abundance distribution of chemicals in 
phase II is close to the equilibrium distribution. Hence, the transition from phase II to I 
is a bifurcation from a branch near equilibrium to that far from equilibrium. It should be 
noted that the power law with respect to the abundance distribution was also observed in 



the previous model, without consideration of energy conversion 12[], while phases II, III, and 
IV were first observed in this study. 
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It is interesting to note that the power law behavior exhibited in the protocell model in 
this study is also commonly observed in the existing cells with regard to gene expression 



levels 1 12[|. In particular it is confirmed in the log phase in bacteria where cells grow efficiently 
as in the phase I in our model. Note that the power law in abundances in protocell models 
is observed for a state in which a cell grows efficiently keeping its chemical composition. 
The distribution of abundances over chemicals may provide a proper statistical measure for 
a cellular state. 

Because of the simplicity of our model, it is not so easy to directly compare our results 
with the behaviors of the existing cells. However, it is still interesting to note the similarity 
of the phases observed in this study with the log, stationary, and dormant phases. The log 
phase appears under a nutrient-rich condition, where the cell has a high growth rate with a 
high influx rate. The stationary phase appears when a nutrient becomes exhausted or toxic 
waste products become excessive, and the cell growth slows. In the dormant phase, the 
cell growth rate is suppressed even under a nutrient-rich condition s|. These behaviors are 
reminiscent of phases I, III, and II, respectively. In this respect, it will be interesting to check 
whether there is hysteresis between the transitions from the log phase to the other phases at 
the single-cell level and to measure the efficiency of cell growth. Mutual inhibition between 
the cell growth and catalyst abundances to support the bistability of the two phases can be 
examined directly. As mentioned above, statistical distribution of chemical abundances is 
distinct by each phase, which may be examined from gene expression analysis. So far, no 
reliable microarray data are available for stationary and dormant phases, which should be 
important in the future. 

The present model does not assume a specific reaction network and intracellular spatial 
organizations. Currently, several attempts are being made to construct an artificial cell that 



grows and replicates |24|, and the basic structures of the present model are often satisfied 
with such an artificial cell. Hence, it will be interesting to search for phases by measuring 
the growth rate and the influx rate for cell growth. 

In general, there is no established theory for thermodynamics far from equilibrium. The 
metabolic efficiency rj discussed above is also not, in a strict sense, the thermodynamic 
efficiency However, by restricting the theory to a steadily growing state, we may hope 
to develop a thermodynamic theory on the basis of which the efficiency of cell growth is 
analyzed 25] . The present model exhibits three distinct phases with regard to the influx 
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and growth rates, with unique metabolic efficiencies. Such a thermodynamic theory, if 
established, will account for these phases and efficiency of cell growth. 
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